The walk-sum method for simulating quantum 
many-body systems 



P-L Giscardi, M Kiffner^ ^ and D Jaksch^ ^ 

^Clarendon Laboratory, Department of Physies, University of Oxford, Parks 
Road, Oxford OXl 3PU, United Kingdom 

^Centre for Quantum Technologies, National University of Singapore, 3 Science 
Drive 2, Singapore 117543 

E-mail: p .giscardlOphysics . ox . ac.uk 



Abstract. We present the method of walk-sum to study the real-time 
dynamics of interacting quantum many-body systems. The walk-sum method 
generates explicit expressions for any desired pieces of an evolution operator U 
independently of any others. The computational cost for evaluating any such 
piece at a fixed order grows polynomially with the number of particles. Walk-sum 
is valid for systems presenting long-range interactions and in any geometry. We 
illustrate the method by means of two physical systems. 
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1. Introduction 

The quality of control over atomic systems in state of the art experiments is such that 
one can now address and observe quantum evolutions of individual atoms in optical- 
lattices [TJ [21 [3l 13] . Additionally, coherent inter-atomic and light-matter interactions 
can be made strong enough to occur on short time-scales compared to incoherent 
processes. This reveals the system's unitary evolution at the individual constituent 
level which is of fundamental interest in the study of quantum many-body phenomena. 
Several applications like quantum simulation and quantum computing schemes also 
rely on this information [5] . On the other hand, the theoretical simulation of such non- 
equilibrium many-body dynamics is very challenging, mainly due to the exponentially 
large number of relevant degrees of freedom. 

This has led to the rise of an active area of research focused on the development of 
theoretical methods to simulate the dynamics of quantum systems. Among the most 
well known techniques are the Density Matrix Renormalization Group method [5], 
related Tensor Network approaches [7] and the Time-Evolving Block Decimation 
method [8]. These techniques have proven to be very successful in their respective 
domains of applicability (e.g. O [101 EDj i-C- mostly ID systems with nearest 
neighbour interactions. On the contrary, and in spite of extensive work in the field, 
very few techniques exist in 2D and 3D [T^ [T3] and for systems exhibiting long-range 
interactions. 
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In this article wc introduce an approach to the real-time evolution of quantum 
systems that is not restricted by the system geometry and works in the presence of 
long-range interactions. Our method, termed walk-sum, describes quantum evolutions 
as resulting from a superposition of all the histories of the system, a discrete 
analog to Feynman path-integrals. We show that with walk-sum one can evaluate 
arbitrarily chosen pieces of the full evolution operator U independently of one another. 
Additionally, we show that each piece can be approximated efficiently, i.e. in a number 
of operations that scales polynomially in the number TV of involved particles. 

This article is organised as follows. In section[2] we introduce the walk-sum 
approach for the simulation of quantum many-body systems. In order to illustrate 
the method and its capabilities, we discuss two physical examples in section|3l The 
computational cost for implementing the walk-sum method is described in sectional 
Finally, we present a conclusion and outlook in section [S] 

2. The walk- sum method 

We consider a quantum many-body system § to be comprised of two parts S and S' . 
The method of walk-sum is based on the following observation: If S' is a large set of 
constituents whose dynamics is frozen, then all interactions between S and S' can be 
evaluated exactly and the evolution of S is described by a small effective Hamiltonian. 
Walk-sum expresses the true many-body dynamics in terms of such simple situations 
with a frozen S' and a few evolving constituents S. This approach might seem similar 
in essence to mean field theory where an atom of interest interacts with a field resulting 
from the mean behaviour of all other particles. The difference is that here we make 
the mapping from many-body to few-body dynamics exact, that is we make S interact 
with all possible fields it could be subjected to depending on the configuration of S' . 

In practice, the walk-sum method divides the matrix representing the time 
evolution operator of a finite-dimensional quantum system into small pieces. These 
pieces are called conditional evolution operators and are introduced in section 12.11 
where we describe the general approach of the method. In section 12.21 we derive an 
equation that allows one to calculate any conditional evolution operator and therefore 
any desired piece of the full evolution operator. We then show in section [231 that the 
generation of conditional evolution operators can be represented by walks on a graph. 
This facilitates the evaluation of the walk-sum expression. In section 12.41 we address 
the evaluation of conditional evolution operators via Laplace transforms. Finally, we 
demonstrate how the method allows one to evaluate expectation values and conditional 
expectation values in section 12.51 

2.1. Model system and general approach 

We consider a quantum many-body system § that is described by a time-independent 
Hamiltonian H and comprised of N physical constituents. The state space of the fcth 
constituent is spanned by dk internal levels, where dk can be different for every k. 
Let {|«fc)} be an orthonormal basis of length dk associated with constituent fc, and 
Pk,i = \ik){ik\ is the projector onto state \ik). Next wc split the system § into two 
parts S and S' . Here S is taken to be constituent s, and all remaining constituents 
are grouped in S". Note that constituent s can be comprised of one or several physical 
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particles. We define so-called projector-lattices 

Acts on S" Acts on S 

e^i^ ®keS'Pk,ik ®1-s= ® Is (1) 

that are projection operators in the state space of the full system §. The index /i 
is a short-hand notation for the configuration onto which S' is projected and Xg is 
the identity operator on S. There is a total of Ofc^s '^fc different projector-lattices 
corresponding to all possible combinations of projectors PkSk ■ Projector-lattices are 
themselves orthogonal projectors and the ensemble of all possible projector-lattices is 
closed, 



(2) 



where I is the identity operator on the state-space of the whole many-body system §. 
We emphasise that our method will work for any set of projector lattices {s^} that 
fulfil Eq. ©. 

Next we partition the time evolution operator U = exp[—iHt] (throughout this 
article we set h= 1) into submatrices with projector lattices, 

£^U(t)e^ = (E)kes'Tk,,k^jk (E> U^^^{t) ^ (g) (3) 

Acts on 5' Acts on S 

Here TkA-i-j — \'>-k)(jk\ denotes a transition operator in subsystem k, and Uu-i^fj,{t) is a 
so-called conditional evolution operator that is represented by a (small) ds x ds matrix. 
Uu^^fi, (t) generates the time evolution of subsystem S provided that (i) the remaining 
many-body system S' was initially in configuration |^), and (ii) the state of S' at time 
t is given by The closure relation in Eq. ([2]) implies that the full time evolution 
operator U(t) can be represented in terms of conditional evolution operators, 

u{t) = J2 { ® two }. (4) 

' ^ Acts on S' Acts on S 

The method of walk-sum allows one to compute conditional evolution operators 
efficiently, i.e. in a number of operations scaling polynomially in the number of involved 
particles, and independently of one another. If in a particular physical problem only 
a small subset of all possible configurations of S' is relevant, then our approach allows 
one to approximate the full time evolution with a moderate computational effort |14j . 



2.2. Derivation of the conditional evolution operators 

In this section we derive an explicit expression for conditional evolution operators that 
allows one to generate them without calculating the full time evolution operator. To 
this end, we write the time evolution operator U(t) of the complete system S as 

M M 

U{t) = cxp[-iHt] = \im Y[SU^ lim H { H ^^}' 

n— 1 n— 1 a 

where SU ~X — iHdt is the infinitesimal evolution operator from (n — l)St to nSt and 
M = t/6t. In the last step in Eq. ([5]) we inserted the closure relation for projector 
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lattices [see Eq. ([2])] . In the following, we derive explicit expressions for the conditional 
time evolution operators Uy^^{t) by a re-summation of terms arising in iyU {t)ifj_. Let 
us first focus on e^U{t)i^, and in particular the term where all intermediate projector 
lattices in Eq. ^ are all equal to i^, 

M 

= n ^^)^^ ^ eJUeJU . . . e^SUe^. (6) 

n = l 

The physical meaning of this term is the following. First S' is projected into 
configuration i/, then 6U evolves the full system § from to St. Then i^, projects 
S' onto configuration v, followed by a time evolution of S for St, and so on. In the 
limit St — > 0, 5" is frozen by continuous (Zeno) measurements of and S evolves 

freely under the influence S' in configuration v. We find 

M 

lim = hm TT (f, - i e.HeJt) = u,{t, 0), (7) 

n=l 

where 

u,{t,t')=e,e-'^-^"'"^'-*'^ (8) 

is a so-called Zeno evolution operator that evolves S between times t and t' while S" 
is frozen in configuration ji^). Zeno-evolution operators obey the following relations 

u^u^ = if v ^ fi, iyUy = Uyiy = Uy, aud iiluy ~ Uyu\, = iy. (9) 

Next we discuss the projection eyU{t)ef^ for [-i^v and focus on those terms that are 
of the form 

= eJUiySU . . . eJUiy SU e^SUe^SU . . . e^SUe^ (10a) 

— —iStiySUiySU. . .iySUiy SyHe^ e^SUif^SU . . .if^SUii^,, (10b) 
' " ' V ' 

= X {M-M' times 5U) =y (M'-i times SU) 

where we employed iySUi^ = —iStiyHe^. Now let t = MSt and t' ~ M'St. The 

term lA^^l can be interpreted as follows. Initially S' is projected into state |^) and 
we observe that lim^(_j.o y = u^(t',0). This is the Zeno evolution operator for S 
that keeps 5" static in state \fi) between and t' . Then the subsystem 5" makes 
a jump from configuration to configuration [term iyHi^], followed by the 
evolution of S between t' and t with S' frozen in state \i/) [term X in Eq. (jlOap fulfils 
Imist-foX = Ufj_{t,t')]. Note that the limit i — >■ cannot be taken directly in ul]} , 
since there are M -> c« ((5t — > 0) similar terms appearing in eyU{t)e^ that differ from 

1^!)]} only in the position of the jump operator iyHi^, i.e. the time t' at which S' 
makes the transition fi ^ i/. Summing over all these terms and taking the limit t — > 
yields 

lim y^-iStX{t,qSt)Hy{q5t-6t,0)^-i Uy{t,t')Huf,{t' ,0)dt' , (11) 
-5*^0^ Jo 

and the integral is thus a continuous sum over the time for the jump fi v to 
occur. The full expansion of £yU{t)i^ also contains terms with an arbitrary number 
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of intermediary eonfigurations 77^ between the initial and final configurations fi and ly, 
respectively. For each of these configurations, the jump rjg — > r/q+i can occur at any 
time tq between and the next jump at time ig+i. Thus 

euU{t)ef,^u^{t,0)S^^^+ (12) 

ft fin ft2 

^i""^ / / ■■•/ U^{t,tn)HUrj^{tn,tn-l)H . . .Urji{t2,ti)Hu^{ti,0)dti . . .dtn, 

where the sum over n counts the number of jumps between different intermediate 
configurations rjq. The result in Eq. (|12p can also be formally derived with graph 
theoretical methods |15j . In addition, we provide a short alternative proof of this result 
in [Appendix A[ demonstrating that Eq. ([T^ can also be regarded as a perturbative 
expansion of the time evolution operator in the transitions underwent by S' . While S 
evolves smoothly in time, Eq. ((T2|) describes 5" as evolving stroboscopically in its state- 
space: e.g. 5" is static in configuration /i from to time ti, then jumps instantaneously 
to another configuration i> where its stays static until t2 , and so on. Each succession of 
configurations adopted by S" from to time t is a possible history of 5". The dynamics 
of S as obtained by Eq. then appears as the superposition of the effects on S of 
all the possible histories of S'. 

The expression for et,t/(t)e^ in Eq. (|I2|) simplifies greatly on using the mixed- 
product property of the tensor product {A ® B){G ® D) = AC ® BD. Indeed, this 
implies that for any two projector-lattices and there exist small ds x ds matrices 
iJj/, H,jj-^ and Uv-i^^it) such that 

Acts on S' Acts on S 

e^Hi^ ^ (g)kes' Pk,r, = ® , (13a) 



®kes' Pk.^,<^^e-'"^'' - \iy){iy\ ® g-^^-*, (13b) 
evHe^= ®keS'Tk.ik^]k®Hi,^^^ ® H^^^,. (13c) 

The Hi, matrices are called statics because they physically correspond to small effective 
Hamiltonians driving S when S' is static in \ v), i.e. frozen in configuration v. Similarly, 
the H^4^^ matrices are called flips and represent how a jump of S' from /i to v affects 
the dynamics of S. The are Hermitian and {H^^f^y = H^^^^. Using these matrices 
and the mixed-product property we can completely separate the evolution of S from 
the evolution of S' . It follows that Eq. ([T2|) is equivalent to 

00 

U,^^{t) = e-'^^*^^,, + (14) 

n=l 

,„ Jo Jo Jo 



X 

r;i,...,r). 



Note that this expression for U^^^{t) only comprises ds x ds matrices. With Eq. (|T4)) 
we achieved to transform the complications associated with the real-time dynamics 
of a quantum many-body system into the sum over the intermediary configurations 
rji, . . .rin of S". In the following section [273l we show that each contribution to the sum 
in Eq. (jl4p can be represented as a walk on a graph. This visualisation clearly brings 
out all the physical processes contributing to U^^^. Most importantly, we show that 
this graphical representation facilitates the evaluation of Eq. substantially. 
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Figure 1: (Colour Online), (a) A square, the graph Q of all configurations available to 
S", here an ensemble of two spin- 1/2 particles labelled 1 and 3. Each edge represents 
a transition allowed by the Hamiltonian in Eq. (|15|) . Highlighted in solid- red and 
dashed-blue are two walks of S" on C/, i.e. two possible histories of S' between time 
and time t. (b) The histories of S' represented by the solid-red and dashed-blue 
walks in (a) where ti, t2, t'l and are the jumping times. This representation of the 
histories can also be thought of as the discrete equivalent of Feynman diagrams, the 
physical processes being the spin flips caused by the Hamiltonian. 



2.3. Representing the time evolution on a graph 

The sum over the intermediary configurations r]q in Eq. (jl4p can be interpreted as a 
sum over the histories of S", a history being a succession of configurations adopted by 
S' from to time t. To visualise the histories, we construct a graph Q as follows: 

i) For each configuration v available to 5", draw a vertex w^. 

ii) For each non-zero flip, i.e. when S' is allowed to make a transition between two 
speciflc conflgurations pi and i^, draw an edge between and v^. 

Now wc can sec that a multiplication by e"*^""^ in Eq. (fT4|) represents the evolution 
of S while S' is stationary at vertex during time r. A multiplication by -ffjy^-^ in 
Eq. (fT4|) corresponds to an instantaneous move of S' along the edge — > Vi, . Conse- 
quently, the histories of S' appear as walks on Q and their superposition is obtained as 
the sum over all the walks, in complete analogy with Feynman path-integrals. If the 
system had continuous degrees of freedom, the sum over histories, that is the sum over 
the walks, would become an integral over histories which corresponds to an integral 
over the walks, i.e. a Feynman path-integral. An equivalent to the Feynman diagrams 
also exists in our situation: this is the succession of physical processes underwent by 
S' during a history. Finally, from a practical point of view, the graph itself facilitates 
the visualisation of the histories of S" contributing to a given conditional evolution op- 
erator. To clarify these notions, we now construct the graph of a small physical system. 

Example. Consider a ID chain of three spin-1/2 particles, and choose S to be the 
central spin, labelled 2. S" is thus comprised of the other two spins 1 and 3. Suppose 
that the Hamiltonian of the system is 

i? = ^Aa:+^JaX+\ (15) 

i=l i=l 

where ax,z are Pauli matrices and A and J are two real parameters. There are 2^ = 4 
orthogonal configurations available to S' and Q has 4 vertices. Indeed, choosing the 
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configurations to be those specifying the direction of tlie spins along the z-axitQ, the 
operator lattices are etita = I tits) (tits I, £tiU = I tits) (tits I, e^its = I tits) (tits I 
and £4.14,3 = I tits) (tits I- We find that the flips and statics are given by 

i/t,t3 =2Ja:+A(7:, i7^,43 = -2Ja: + Aa^, (16a) 
Huts = H^iU = H,^^^ = AIs, (16b) 

where v and /x are two configurations of S' that differ by exactly one spin flip. The 
graph Q has therefore 4 edges and is the square. The graph Q as well as two examples 
for histories of 5', two walks and the equivalents of Feynman diagrams are all repre- 
sented in Fig. [T] 

In order to make the superposition of histories manifest in the expression for 
U^i-fi(t), we relabel the sum over intermediary configurations into a sum over all the 
walks on Q. If the initial and final configurations of S' appearing in Eq. ((T4|) are 
respectively fi and h', then the relevant walks all start on vertex Vf^ and terminate on 
the vertex Vi,. The conditional evolution operator C/j/<_^ is therefore 

00 

Jo Jq Jo 

where Wg-,^^-n is the ensemble of walks of length n on ^ from to v^. Eq. (|17p is the 
central result of this work. The sum over the walks clearly identifies the successions of 
configurations of 5" that are allowed by the Hamiltonian. Finally, Eq. ([TT]) also presents 
a mathematical advantage. It points to the general approach to matrix functions from 
which the walk-sum method ultimately stems. This approach relies on graph-theoretic 
arguments with walks on graphs playing a central role |15] . 

2.4- Evaluation of conditional evolution operators 

Conditional evolution operators are most conveniently evaluated in the the Laplace 
domain, because this transformation turns the nested convolutions in Eq. (|17p into a 
product of matrices. Let M^{s) = £[exp(— ii/^i)] = {sis + iH^)^^ be the element- 
wise Laplace transform of exp(— iiJ^t) with s the Laplace domain variable. It follows 
that the Laplace transform of a conditional evolution operator is given by 

00 

Uu^^is) = M,{s)5^^, + Y. ^^-(^) • ■ ■ ^■"1^'^ ^^^(*)- (18) 

This is the expression we use in computations. It only involves multiplications and 
additions of da x ds matrices. Considering the case where all the M^(s) can be obtained 
analytically, we remark that the matrix elements of the Mfj_{s) are ratios of polynomials 
in s. The roots of the denominator polynomial are — iA^ with an eigenvalue of 77^. 
This is therefore also true for any element of the sum in Eq. ()18p and performing the 



I One could equally well choose to specify the spins along another axis. 
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inverse transform into the time-domain ean be done aceording to the formula 

C/.^mW = ^Res[e^*;7,^^(s),-zA^], (19) 

where Res[/(s),A] denotes the residuum of /(s) at A. In aU practical situations the 
infinite sum in Eq. (jl7p is truncated at some order K. With this approximation 
only walks of length n < K are taken into account. This approximation is justified 
because the series in Eq. ()17|) is absolutely and uniformly convergent for any finite time 
provided § comprises a finite number of particles. In order establish this result, we 
consider the contribution ul^fj,{t) of all walks whose length is exactly K to Uyi-^(t). 
A bound for the norm of this quantity is given by 

\\Uf^,m<<^N{K)^^^^Q (K^^), (20) 

where = max„{||i7,,^^j.<_,,^ ||}, || • || is the spectral norm and a^iK) is the number of 
walks of length K on the graph Q. In |Appendix C| we show that a^iK) is a polynomial 
in the number of particles N of the system, and hence ||?7i^^(t)|| tends to zero as 
K ^ oo. With the help of the inequality in Eq. ([^0]) it is straightforward to show 
that the sequence Uk = ^n=Q u!^fj.{t) is Cauchy, which establishes the convergence 
of the walk-sum series in Eq. (|17l) . This result follows also from the convergence of 
the Dyson series for finite-dimensional matrices |16j and its relation to the walk-sum 
method (see [Appendix A[ ). 

Several conclusions can be drawn from Eq. (pUj) . First, it shows that the evaluation 
of walk-sum series is free from the numerical sign problem. Second, an increase in 
time requires to include higher order terms in Eq. (jl7p . Third, the convergence in the 
evaluation of U^-i^n can be established if an increase of the order K does not alter the 
result within a given accuracy. All these conclusions are confirmed by our experience 
with example systems in section |3] and |14j . 

A more detailed analysis of the convergence of walk-sum is beyond the scope 
of this paper. However, we emphasise that we have always observed that walk-sum 
converges much faster with the order K than Eq. (|20| suggests. Finally, we note that 
conditional expectation values converge much faster than ordinary expectation values 
because the former quantities are insensitive to the norm of the overall state vector. In 
the next section we discuss how expectation values and conditional expectation values 
can be calculated via conditional evolution operators. 

2.5. Expectation values and conditional expectation values 

Here we demonstrate how conditional evolution operators can be employed for the 
evaluation of expectation values and conditional expectation values. For simplicity 
we assume that the initial state of the full system is given by |V'(0)) = Im) ® IV's(O)), 
where |'!/'s(0)) is the initial state of ^f^. The state vector of the full system can then 
be obtained via Eq. (|4]) and is given by 

\m) = u{t)\m) = E {i^) ® mt)).^^], (21) 

§ All the relations of this sections are easily extended to more general initial pure states ® 
|i/'s{0)) upon introducing a sum over ^l as appropriate. 
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where \ips{t))u^fi ~ t^i/^/j(i)|V's(0)} is a piece of the wavefunction of the fuh system. 
Any mean value of an observable O can be caleulated via Eq. (PT|) . With the definition 
e^'Oii, = ® Oiy'i,, we find 

(O) = {Am {Uu'^t^y O.v U,^^ IV.(O)). (22) 

It follows that any expectation value of an observable O can be computed directly 
from conditional evolution operators if O is expanded in the projector lattice basis e^. 

Next we investigate in more detail the meaning of the piece \'>ps{t))^4-fj, = 
U„^^{t)\'ips{0)) of the full state vector. To this end, we write the projection of 
onto ii, as 

\^it)) = \A{t))u^^. = £ WQ^ (23) 

£=1 

where the probability amplitudes ae of state \iyis) (S" in configuration v, S in state \is)) 
are determined by U,^j^i^i,{t). This shows that U^^^{t) directly evolves the projection 
of the full state vector \^{t)) onto a dimensional subspace corresponding to 5" in 
configuration v. We emphasise that conditional evolution operators are submatrices 
of U [see Ec^. (jl])] and therefore not necessarily unitary. It follows that Uu^^i{t) does 
generally not conserve the norm of 1-08(0)). Indeed, Eq. (|23|) implies that the norm 
of |'0s(i))iy<_p is given by the probability of finding S" in configuration v at time 
knowing that it was initially in configuration /i, 

(0,(0)1 iU^^^it))^ U.^^it) = {e.)t- (24) 

We now turn to conditional expectation values and show that some of them can be 
computed directly via a single conditional evolution operator. To this end, we note 
that the expectation value of an arbitrary operator Os acting on S with respect to 
Hs{t))v^t, is given by 

(0,(0)1 {U,^^,{t))^ Os U,^^,{t) 10^,(0)) = (0.e,)t, (25) 

where {Os£v)t is the expectation value of Osiy at time t. Eqs. ([24l[25|) directly yield the 
conditional expectation value (OsM of Og, knowing that S' , initially in configuration 
fi, is in configuration v at time t, |j 

_ {Osi.) _ (0s(o)|(^.^^(O)^o.t/.^^ (010.(0)} ^2^^ 

(0;,(O)|(C/,^^(t))V,^^(t)|0;,(O)) 

Finally, Eqs. ([2]) and (|26|) imply that the full set of conditional probabilities allows one 
to compute the mean value of an operator acting on S, J2,yi^s/'^){£iy) = (Os)- We 
illustrate the preceding results in section [31 where we apply the walk-sum method to 
two physical systems. 



3. Examples 

In this section we present two examples for the walk-sum method. First we consider a 
simple case where 3 spins evolve according to a nearest neighbour interaction and 

II The fact that (Osii,) / (ii,) is indeed the conditional expectation value follows from [Os,ei/] = 
and £^ = ii^ | 17| . 
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Figure 2: (Colour online). Top figure: evolution of (t2 / tits) evaluated by walk- 
sum at order 2 (solid blue line) compared to the exact solution (dashed red line). 
Bottom figure: evolution of (t2 / tits) evaluated by walk-sum at order 4 (solid black 
line) compared to the exact solution (dashed red line). Note that order 4 is the next 
non-zero order after order 2. Parameters: J/A — 4. 

we compare the walk-sum predictions to the exact results. Second, we consider 
the dynamics of a ID chain of spin-1/2 particles interacting via strong long-range 
interactions. 




1 2 3 4 5 6 7 

Jt/n 



3.1. Dynamics of 3 spins 

Consider a small ID chain with 3 spin-1/2 particles evolving according to the 
Hamiltonian in Eq. ()15p . Suppose that initially all three spins are up along z, i.e. 
IV'(O)) ~ I tit2t3) a-nd suppose that we are interested in the time evolution of the 
conditional expectation value (t2 / tits)- This is the probability for spin 2 to be 
up along z knowing the other 2 spins are up along z. To implement walk-sum, we 
take S to be the spin 2 and S' comprises spins 1 and 3. Similarly to the example of 
section 12. 3| we choose the configurations of S' to specify the direction of these spins 
along the z-axis and obtain the flips and statics of Eq. (|16p . The graph Q is the square 
represented in Fig. [TJ According to Eq. (pS)) . the quantity (t2 / tits) is given by 

i^s (0) I (Uu^.^m.)^ Pu Uuu^u^. I (0)) 

\\2 I tits) - ; -f , 

where P-^^ — \ t2)(t2 | and \iljsiQ)) = | t2) is the initial state of subsystem S. 
Considering the closed walks off vertex f-t-j-i-, on Q and keeping only walks of length 2 
or less, the conditional evolution operator in the Laplace domain U^-^^^^^^-^^^ is given 

by 

Walk of length ^wo walks of length 2 



«^t2t3<-t2T3(s) - A^Tit3 + (-^A)2 Af^,t3(A/;,t3+Mtii3)A'%T3, (28) 

where the factor (— lA)^ comes from the flips. The two walks of length two taken 
into account here are those highlighted by solid-red and dashcd-blue lines in Fig. [1] In 
Fig.[2]we compare the walk-sum result for the time evolution of (t2 / tits) to the exact 
result. For short times Jt < 37r, the agreement between the walk-sum estimation of 
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(t2 / tits) at order 2 and the exact result is excellent. Differences appear for Jt > Stt 
and grow larger with time. At order 4, which is the next non-zero order, the difference 
between the walk-sum estimation and the exact result does not exceed 5 x 10"'^ for t 
up to Jt ~ Ttt, which represents a relative error of less than 0.5%. 



3.2. Dynamics of a ID spin chain with long-range interactions 

In our second example we demonstrate that the walk-sum method allows one to 
investigate the quantum dynamics of systems with long-range interactions. More 
specifically, we consider a one-dimensional chain comprised of 27 spin-1/2 particles 
that is described by the Hamiltonian 

{AP-il)< + iE^^^^1' (29) 

where A, Q and U arc real parameters and we introduced the Pauli matrix cr* of 
spin i. Throughout this section, we choose the z axis as the quantisation axis and 
take the eigenstates | ti) and | ii) of the Pauli matrix <tI as basis states of spin i. 
Furthermore, we set P' = |ti)(t» \ ^ + c^i)/2 and G' = \U){U \ = ~ 
The term proportional to U in Eq. (|29p describes a long-range interaction between the 
spins that scales as XjP? with the distance R between two spins. 

We now study the real-time dynamics of the conditional expectation value 
(tc/tj) = (tctj)/(tj), where (tctj) ((tj)) is the probability for spins c and j (spin f) 
pointing upwards, (tc/tj) describes the probability that spin c at the centre of the 
lattice is in state |tc), knowing that spin j ^ c is pointing upwards as well. For the 
implementation of walk-sum, we choose S' = c to be the spin at the centre of the lattice 
and assume that all spins are initially prepared in \\). We choose the projector-lattices 
to be tensor products of projectors P* and G* on S", 

e\...»„ = «)res'\{»i...»„}G^ P"-' ®Is- (30) 

The projector lattice Si-^,,,i„ projects S' in a state where all spins at positions ii . . . i„ 
are pointing upwards, and all other spins of S' are pointing downwards. With this 
definition, we find that the flips are given by 

Hii...i„±i-^ii...i„ = ^2 ^-^s' and otherwise, (31) 

that is H^^^ 7^ if and only if v differs from fi by exactly one spin flip. The evaluation 
of the statics yields t 

(ip,,f)e{ii.-.i„} 

The conditional expectation value (tc / tj) = (tcti)/(ti) is composed of the 
probabilities (tctj) and {tj)- These quantities can be computed via Eq. ([22]) . 

itsW = I E iUvU)^uyP-' C/,(,)^iJ;.), (33a) 
vU) 

iW = iU I E iUnU)^uy UnU)^u\is). (33b) 
vU) 

In this equation, the sum runs over all configurations ri{j) of S' where spin j is pointing 
upwards. We compute the conditional evolution operators in the Laplace 
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Figure 3: (Colour online). The conditional probability (ts/tj) after a 7r-pulse fit ~ tt 
with A = and U = —20051. Subsystem S is represented by the central spin (label 
0). Insets show the time evolution of (ts/tj) for j = —13, — 2, 4 and 6. 



domain via Eq. (|18p . For example, the conditional evolution operator Uji-o with only 
spin j pointing upwards in the final configuration 77 (j) of S" and for walks of length 
n < 3 is given by 

;7,^o(s)^yA7,Mo+ (^y) ^2 M,(Mo + M,,p)(Mp + M,)Mo. (34) 

In this expression, the factors in/2 and {irt/2)^ come from the flips and Mij...i^(.s) = 
£[exp(— i_ffij...i^i)] , where the static Hi^,,,i^ is deflned in Eq. (|32]) . The sum over 
p in Eq. p4[) accounts for an up-spin in an intermediary configuration of 5" that 
is not present in the initial and final configurations. Note that intermediary spin- 
excitations are similar to virtual particles appearing in the Feynman diagrams of QED. 
This observation further substantiates the analogy between the walk-sum method and 
path-integrals. 

In the evaluation of Eq. p3p we consider all configurations 77 (j) of 5" with at 
most 4 up-spins. For each conditional evolution operator we took into account all 
walks of length 4 or shorter, i,e., we truncated the series in Eq. ([T8| at n = 4. We 
find that these approximations are justified because the result for (t,s/tj) converges 
with the terms taken into account and for the chosen parameters. In particular, we 
consider the regime [/ ^ 51 where the interaction U between nearby spins dominates 
the coherent driving term 51crj/2 in Eq. ([29|) . In this case, the strong interaction 
inhibits the simultaneous excitation of two nearby spins, and only configurations of S' 
with a few, evenly distributed up-spins will contribute significantly to Eq. ([55)1 . 

The final result for (ts/tj) in the time domain is shown in Fig. ([3]). Note that 
(ts/tj) is much smaller than unity for \j\ < 6. It follows that the excitation probability 
for the central spin at site (subsystem S) is very small if another spin at site \i\ < 6 
is pointing upwards. This result is characteristic of a blockade effect [HI [TH [20] that 
is well known for Rydberg atoms coupled by the dipole-dipole interaction. Indeed, we 
show in [Appendix B| that a physical realization of the Hamiltonian in Eq. ([^H)) can 
be achieved with a one-dimensional chain of atoms in their ground states that are 
subsequently laser-excited to Rydberg states. 

4. Computational cost: polynomial scaling in A'^ 

In this section we show that the number of fioating point operations required to 
approximate a conditional evolution operator to order K scales polynomially in the 
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number of particles N of the system. For simplicity and without loss of generality, 
we assume that all particles have d internal levels and that S contains p particles, 
and thus ds = d^. The evaluation of a single walk of length if to a conditional 
evolution operator requires 2K + 1 multiplications of ds x ds matrices (K statics 
and if + 1 flips), i.e. ~ "^Kdg operations. Then all the walk contributions must 
be added together. This represents an additional \Wg-„f_i-K\dl operations, 
being the number of walks of length K on Q between vertices Vi, and u^. We show 
in [Appendix C| that \Wg-^fj,-K\ < poly(A'^) is always bounded by a polynomial in N. 
This is a direct consequence of the sparsity of quantum Hamiltonians: there are not 
so many configurations accessible to S' through K jumps. An upper bound for the 
number of operations Mk required to obtain C/^<-y at order K is therefore given by 

J\fK < {'iKdl + dl) poly(iV). (35) 

It follows that the computational cost for the evaluation up to order K of a single 
conditional evolution operator scales polynomially with the number of particles TV. 
On the contrary, the progression of the accuracy with the order K and the number of 
particles N are open questions. However, we find that the accuracy of walk-sum at a 
given order K exhibits a very weak dependence on N in the case of the example in 
section 1321 This observation suggests that the walk-sum method scales polynomially 
in the number of particles N for strongly interacting spin chains where blockade effects 
limit the total number of excitations in the lattice. 

5. Conclusion 

We have introduced the walk-sum that allows one to investigate the real-time dynamics 
of quantum systems. We have shown that walk-sum can approximate any desired 
piece of an evolution operator independently of any other. Furthermore, the number 
of operations involved in approximating any such piece at order K scales polynomially 
with the system size. This result holds independently of the system geometry and the 
nature of its interactions. We find that the walk-sum method is very well suited for 
the study of strongly interacting spin chains where blockade effects limit the number 
of final excitations in the lattice. Note that a physical realisation of these systems is 
given by Rydberg atoms interacting via the long-range dipolc-dipole interaction. These 
systems are well confined to a part of their state space whose size grows polynomially 
with the number of particles. In addition to the efficient approximation of single 
conditional evolution operators, these systems allow one to approximate the entire 
evolution operator with moderate computational effort. The reason is that the number 
of all conditional evolution operators contributing substantially to the norm of the 
state vector scales polynomially with the number of particles N . 

Walk-sum is similar to Feynman path-integrals in that it describes a quantum 
evolution as resulting from the superposition of all the possible histories of the system. 
It differs from path-integrals in that it is well adapted to systems with discrete degrees 
of freedom and can independently evolve any chosen sub-ensemble S of an otherwise 
large many-body system. We will show elsewhere that walk-sum can be extended 
to systems with continuous degrees of freedom. In this case, and if S has only one 
internal level, then the walk-sum method is formally equivalent to path-integrals. 

Walk-sum is also well suited to describe the effect of post-selection in quantum 
systems. Indeed, since a conditional evolution operator specifies both the initial and 
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final configurations of S", walk-sum is really a time-symmetric description of quantum 
mechanics. In particular, we will show that the Aharonov-Bergmann-Lebowitz rule 
of the Two State Vector Formalism [5T] follows from a simple constraint on the 
walks of the graph Q and is straightforwardly demonstrated with conditional evolution 
operators. 

Finally, it is worth noting that the deep origin of the computational speed- 
up achieved by walk-sum resides in exact summations of infinite families of terms 
performed at the graph-theoretic level by Eq. (flTl [T8| . These resummations make 
walk-sum fundamentally different from power series and are explicitly apparent in the 
full graph theoretic treatment provided in |15| . In fact, further resummations are 
possible that reduce the sum over walks into a sum over simple paths, a simple path 
being forbidden to visit any vertex more than once. The resulting method, called 
path-sum is nonperturbative and computes any conditional evolution operator exactly 
in a finite number of operations. Path-sum can be also seen as a generalisation of 
the method of resolvents. Its use for the study of quantum dynamics will be exposed 
elsewhere. 
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Appendix A. Alternative derivation of Eq. (|12p 

Here we provide an alternative derivation of Eq. (fT2)) for £yU{t)£^ via the well-known 
expression (see, e.g., Complement AI.2 in [22] ) 

oo 

U{t, 0) = exp[-i(i7o + V)t] = exp[-iHot] + ^ (A.l) 

n=l 

/""../" 'e-'-f^"(*-*"Ve-^-f^°(*"-*"-iV...e-*^°'idti...dt„, 
JO Jo 

which is used in time-dependent perturbation theory. In order to establish Eq. (|T2|) 
via Eq. (|A.ip , we expand the Hamiltonian H in terms of the complete set of projector- 
lattices in Eq. © and obtain H ~ Ho + V, where 

Ho^J2 ^-^^-^ ^^Y. ^-^^t^- (^-2) 

If these expressions are plugged in Eq. (jA.ip , the orthogonality of the projector lattices 
directly yields the result in Eq. (fT2)) . This indicates another interpretation for the result 
of Eq. (fT2l) : the transitions undergone by S' can be seen as perturbations affecting the 
evolution of S. 

Appendix B. Rydberg-excited Mott insulators 

In this section we discuss a physical realization of the Hamiltonian in Eq. ([25]). We 
consider a one-dimensional chain of atoms that are prepared in the Mott-insulating 
phase with unit filling via a strong optical lattice. An external laser field couples 
the ground state \g) of each atom to a Rydberg state |r). We assume that atomic 
transitions to other internal levels are negligible such that the atoms are well described 
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by two-level systems. Two Rydberg atoms at lattice site i and j interact via the long- 
range dipole-dipole interaction that is of the form [53] 

A, = (iTTeoRlr' [t^i.fij - 3R;^\fii.Rij)ifij.Rij)] , (B.l) 

where Rij is the relative distance between atoms i and j, and fii is the dipole 
moment of atom i. We assume that all dipole moments are identical and oriented 
along the axis of the lattice. With = /.t^ = llAtJI, the dipole-dipole interaction 
can be written as Aij = — 2/x^(47reoi^)~^/|i — where L is the lattice constant. 
It follows that the chain of laser-driven atoms is described by the Hamiltonian in 
Eq. ((5^ if we set U = — 2/i^(47reoi'^)~^ and if we employ the frozen gas and rotating- 
wave approximations. In this case, A describes the laser detuning with the \g) |r) 
transition, is the Rabi frequency of the driving laser field, and the spin states | t) 
and I D correspond to \r) and \g), respectively. The quantity (ts/tj) describes the 
probability for atom s being in its Rydberg state |rs) knowing that atom j is excited 
as well. In particular. Fig. ([3]) corresponds to the following physical parameters: time 
Ut = TT, lattice spacing L = 4/im, laser detuning A = 0, Rabi frequency = 5MIIz 
and dipole moment fi ~ 2350eao (Rydberg level principal number w 40). 

Appendix C. An upper bound for the number of walks 

Here we derive an upper bound for \ Wg-i,fj,-n\ the number of walks of length n between 
two vertices and Vi, on a graph Q. We employ the same notation and definitions as 
in section m Suppose that the flips of the Hamiltonian 7J of § [see Eq. ((T3)) ] are only 
non-zero for those configurations fi and v of S' where at most q particles change their 
state in a transition from fi to v. For example, the expressions for the flips in Eqs. (|16p 
and (pij) demonstrate that the Hamiltonians in Eqs. (|15p and ((23) have q = 1. In fact, 
most of the common Hamiltonians have q = 1, 2. With p particles in S there are N —p 
particles in S", and thus there are (J^^p^ ways to choose which q particles among N —p 
undergo a transition. Furthermore, for each particle that changes a state there are at 
most d—1 possible basis states available for the transition. In total, there are thus at 
most V < {^q^){d — 1)'' configurations directly accessible to S" from any given one. 
By construction of G, this is also an upper-bound on the number of vertices that share 
an edge with any given vertex on Q. It follows that there are at most 1^" walks of 
length n attached to any vertex of G, and thus we find 

|W^g;.^;„|<(V)"(^-l)'"' 

which scales polynomially in N—p for small q. Note that this upper-bound is very loose 
and we only use it to demonstrate that \ Wg-^i^i-n \ is always bounded by a polynomial in 
N. For example, for a system with + 1 spin-1/2 particles, one particle in S* (p = 1) 
and a Hamiltonian with q = 1, the graph Q is the A^-hypercube. In this case, we find 
that the number of walks of length n = £-\-2k>0, k^N,0<£<N, between two 
vertices separated by a distance £ is 

where pn = N mod 2 and [-J is the fioor function. This is substantially smaller than 
the iV^+2'= of the bound Eq. (|aT|) . indeed for iV > 1, > n, aAr(£ 2k) ~ 2-''{N - 
£/2-lf{£ + 2k)\/k\ < N^+^'' and for n > A^, aN{£ + 2k) ~ 2"^+'^ N^+'^^ < 7v«+2fc^ 
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